

%% load theta
load('../../baseline_model/parameter_estimates_round2_replication.mat')
[~,use] = min(l0);
thetaStar = theta(:,use); %#ok<NASGU>
clearvars -except thetaStar specif


%% load original sample for normalization
load('../../baseline_model/data.mat')
m_gdp=mean(data.gdp_pc);
sd_gdp=std(data.gdp_pc);
m_pop=mean(data.pop);
sd_pop=std(data.pop);
m_polity=mean(data.polity);
sd_polity=std(data.polity);
m_terrain=mean(data.terrain);
sd_terrain=std(data.terrain);
dd=[data.USdist,data.UKdist,data.FRdist,data.RUdist,data.CHdist];
m_dist=mean(dd(:));
sd_dist=std(dd(:));


%% load everything of interest
load('../../baseline_model/sim_eq_round2.mat')
load('data_decade.mat')
U0=repmat(U0,1,1,5);
intEq=cell(NS/5,N*5);
VR=zeros(size(V_R,1),NS/5,N*5);
VUS=VR;
VUK=VR;
VFrance=VR;
VRussia=VR;
VChina=VR;
for t=1:5
    intEq(:,N*(t-1)+(1:N))=intervEq(NS/5*(t-1)+(1:NS/5),:);
    VR(:,:,N*(t-1)+(1:N))=V_R(:,NS/5*(t-1)+(1:NS/5),:);
    VUS(:,:,N*(t-1)+(1:N))=V_US(:,NS/5*(t-1)+(1:NS/5),:);
    VUK(:,:,N*(t-1)+(1:N))=V_UK(:,NS/5*(t-1)+(1:NS/5),:);
    VFrance(:,:,N*(t-1)+(1:N))=V_France(:,NS/5*(t-1)+(1:NS/5),:);
    VRussia(:,:,N*(t-1)+(1:N))=V_Russia(:,NS/5*(t-1)+(1:NS/5),:);
    VChina(:,:,N*(t-1)+(1:N))=V_China(:,NS/5*(t-1)+(1:NS/5),:);
end
U0=U0(:,:,nmissing);
intervEq=intEq(:,nmissing);
V_R=VR(:,:,nmissing);
V_US=VUS(:,:,nmissing);
V_UK=VUK(:,:,nmissing);
V_France=VFrance(:,:,nmissing);
V_Russia=VRussia(:,:,nmissing);
V_China=VChina(:,:,nmissing);
clear t intEq VR VUS VUK VFrance VRussia VChina
[DMES,intervEq] = createDMES(intervEq,V_R,V_US,V_UK,V_France,V_Russia,V_China,YY,specif);


%% normalize continuous variables
data.gdp_pc = (data.gdp_pc-m_gdp)/sd_gdp;
data.pop = (data.pop-m_pop)/sd_pop;
data.polity = (data.polity-m_polity)/sd_polity;
data.terrain = (data.terrain-m_terrain)/sd_terrain;
dd = [data.USdist,data.UKdist,data.FRdist,data.RUdist,data.CHdist];
dd = (dd-m_dist)/sd_dist;
data.USdist = dd(:,1);
data.UKdist = dd(:,2);
data.FRdist = dd(:,3);
data.RUdist = dd(:,4);
data.CHdist = dd(:,5);
clear dd m_gdp sd_gdp m_pop sd_pop m_polity sd_polity m_terrain sd_terrain m_dist sd_dist


%% estimates of mean utilities
X=[ones(sum(nmissing),1) data.terrain(nmissing) data.gdp_pc(nmissing) data.polity(nmissing) data.pop(nmissing)]; K=size(X,2);
Z_US=[data.USdist(nmissing) data.USally(nmissing) data.UScolony(nmissing) data.USwar(nmissing)]; KZ=size(Z_US,2);
Z_UK=[data.UKdist(nmissing) data.UKally(nmissing) data.UKcolony(nmissing) data.UKwar(nmissing)];
Z_France=[data.FRdist(nmissing) data.FRally(nmissing) data.FRcolony(nmissing) data.FRwar(nmissing)];
Z_Russia=[data.RUdist(nmissing) data.RUally(nmissing) data.RUcolony(nmissing) data.RUwar(nmissing)];
Z_China=[data.CHdist(nmissing) data.CHally(nmissing) data.CHcolony(nmissing) data.CHwar(nmissing)];
Z=[Z_US,Z_UK,Z_France,Z_Russia,Z_China];
reb={1:K,2:KZ}; mpow={[1,3:K],1:KZ};
Y=arrayfun(@(n)find(sum(YY==table2array(outcomes(n,5:end)),2)==size(YY,2)),find(nmissing));

options=optimset('Display','off','HessFcn',@(x,lambda)HessPhat(x,lambda,Y,X,Z,YY,V_R,V_US,V_UK,V_France,V_Russia,V_China,U0,intervEq,DMES,reb,mpow));

% "long" replication
% load('decade_replication.mat','SV','srngSV')
% rng(srngSV)
% theta0=thetaStar+[zeros(length(thetaStar),1),randn(length(thetaStar),SV)/10]; % starting values
% theta=theta0; 
% l0=1:SV; % maximized likelihood values across starting values
% exitn=1:SV; % Knitro exit flags across starting values
% f0 = exitn;
% time=0;
% parfor sv=1:SV
%     tic
%     f0(sv) = Phat(theta0(:,sv),Y,X,Z,YY,V_R,V_US,V_UK,V_France,V_Russia,V_China,U0,intervEq,DMES,reb,mpow);
%     if ~isfinite(f0(sv))
%         theta(:,sv) = NaN(length(theta0(:,sv)),1);
%         exitn(sv) = NaN;
%         l0(sv) = NaN; 
%         time=time+toc;
%     else
%         [theta(:,sv),l0(sv),exitn(sv)]=knitromatlab(@(x)Phat(x,Y,X,Z,YY,V_R,V_US,V_UK,V_France,V_Russia,V_China,U0,intervEq,DMES,reb,mpow),...
%             theta0(:,sv),[],[],[],[],[],[],[],[],options,'knitro1.opt');
%         time=time+toc;
%     end
% end

% "short" replication
load('decade_replication.mat')
[~,use] = min(l0);
time=0;
[theta(:,use),l0(use),exitn(use)]=knitromatlab(@(x)Phat(x,Y,X,Z,YY,V_R,V_US,V_UK,V_France,V_Russia,V_China,U0,intervEq,DMES,reb,mpow),...
    theta(:,use),[],[],[],[],[],[],[],[],options,'knitro1.opt');
time=time+toc;


%%
clearvars -except exitn K l0 N SV theta theta0 time X Y YY Z reb mpow srngSV specif nmissing
save('decade.mat')


%% verify optimality
[~,use] = min(l0);
thetaStar = theta(:,use);
load('decade_replication.mat')
[~,use] = min(l0);
thetaStar0 = theta(:,use);
disp(' ')
disp('verify country-decade results:')
disp(['max. rel. diff. between estimates = ',num2str(max(abs(thetaStar0-thetaStar)./abs(thetaStar0)))])



